function [ mu, Sigma ] = ep( K, pref, gamma )
m = size(pref,1);
n = size(K, 1);

Sigma = eye(n,n);
mu = zeros(n,1);
one_1 = [1;-1];
one_2 = [1,-1;-1,1];

k_inv = inv(K);

converged = false;
iteration = 1;
while (~converged && iteration < 10)
    iteration = iteration + 1;
    sig_prev = Sigma;
    mu_prev = mu;
    
    sig_tilde = zeros(2,2);
    mu_tilde = zeros(2,1);
    
    mu_prev = mu;
    sig_prev = Sigma;
    
    i=1;
    for i = 1 : m 
        sig_t = [Sigma(pref(i, 1),pref(i, 1)), Sigma(pref(i, 1),pref(i, 2));Sigma(pref(i, 2),pref(i, 1)),Sigma(pref(i, 2),pref(i, 2))]
        mu_t = [mu(pref(i, 1)); mu(pref(i, 2)) ]
        
        if(det(sig_tilde) > 0)
            inv_sig_tild = inv(sig_tilde);
        else
            inv_sig_tild = zeros(2,2);
        end
        
        sig_sub = inv(inv(sig_t) - inv_sig_tild);        
        mu_sub = sig_sub * (sig_t \ mu_t - inv_sig_tild * mu_tilde);
        
        sig_star = (gamma^2  + sig_sub(1,1) + sig_sub(2,2) - sig_sub(1,2) - sig_sub(2,1));
        s = (mu(pref(i, 1)) - mu(pref(i,2))) / sig_star;
        coef = (normpdf(s) / normcdf(s));
        if(normcdf(s) == 0)
            coef = abs(s);
        end
        temp = (coef/sig_star) * sig_sub * one_1;
        sig_hat = sig_sub - sig_sub * (temp * temp' + s * (coef/(sig_star^2) * one_2)) * sig_sub;
        mu_hat = mu_sub +  temp;
        
        sig_tilde = (inv(sig_hat) - inv(sig_sub));
        while(det(sig_tilde) <= 0)
            sig_tilde = sig_tilde + eye(size(sig_tilde));
        end
        sig_tilde = inv(sig_tilde);
        mu_tilde = sig_tilde * (sig_hat \ mu_hat - sig_sub \ mu_sub);        
    
        sig_bar = zeros(n,n);
        sig_bar(pref(i, 1),pref(i, 1)) = sig_tilde(1,1);
        sig_bar(pref(i, 1),pref(i, 2)) = sig_tilde(1,2);
        sig_bar(pref(i, 2),pref(i, 1)) = sig_tilde(2,1);
        sig_bar(pref(i, 2),pref(i, 2)) = sig_tilde(2,2);
        
        mu_bar = zeros(n,1);
        mu_bar(pref(i, 1)) = mu_tilde(1);
        mu_bar(pref(i, 2)) = mu_tilde(2);

        while(det(sig_bar) <= 0)
            sig_bar = sig_bar + eye(size(sig_bar));
        end
        sig_bar = inv(sig_bar);
        Sigma = inv(k_inv + sig_bar);
        mu = Sigma * sig_bar * mu_bar;
    end
    
    fprintf('iteration: %d' , iteration);
    norm(Sigma-sig_prev)
    norm(mu - mu_prev)
    
    if (norm(Sigma-sig_prev) < .1 && norm(mu - mu_prev) < .1)
        converged = true;
    end
end
end

